% estimate impulse responses to government spending shock 
% recursive identification; sensitivity analysis
% this version: may 2011

clear all; clc; close all;

load us_quarterly       % generated by getdata.m

saveirf      = 1;

rand('state',1);
nptsvar     = 30;       % number of points of IRF in VAR
teststatio  = 0;        % tests wether VAR is explosive
confidence  = 1.645;    % (1.645: 90 percent); % or 1.96 %(95 percent);
ndraws      = 100;      % number of replications in bootstrap
nlagq       = 4;        % number of lags
beg_samp    = (82-47)*4+1; 
end_samp    = (107-47)*4+4;

disp(['Sample: ' num2str(sampleq(beg_samp+nlagq)) ' - '  num2str(sampleq(end_samp)) ' (dependent variables)  '  ])
gy = mean(g(beg_samp:end_samp)./y(beg_samp:end_samp));
cy = mean(c(beg_samp:end_samp)./y(beg_samp:end_samp));
iy = mean(inv(beg_samp:end_samp)./y(beg_samp:end_samp));

for sensi =1:4;

    if sensi == 4
        beg_samp    = (74-47)*4+1; 
    end
    for jrun = 1:3;

        if jrun==1;
            namshort = char('g','y','c','lr','rx','dy','inf');       
            xdata = [log(g) log(y) log(c) lrr log(rx) dy  inflation];   
            if sensi==3 || sensi == 4;
                xdata = [log(g) log(y) log(c) tbill log(rx) dy  inflation];   
            end
            yshare = [gy 1 cy 1 1 1 1];
        elseif jrun ==2;
            namshort = char('g','y','i','lr','rx','dy','inf');       
            xdata = [log(g) log(y) log(inv)  lrr log(rx) dy inflation];
            if sensi==3 || sensi == 4;
                xdata = [log(g) log(y) log(inv)  tbill log(rx) dy inflation];
            end
            yshare = [gy 1 iy 1 1 1 1];
        elseif jrun ==3;
            namshort = char('g','y','nx','lr','rx','dy','inf');      
            xdata = [log(g) log(y) nx  lrr log(rx) dy inflation];
            if sensi==3 || sensi == 4;
                xdata = [log(g) log(y) nx tbill log(rx) dy inflation];
            end
            yshare = [gy 1 1 1 1 1 1];
        end

        %% specify deterministics
        [obs nvar]=size(xdata);
        exo = []; 
        for ik=1:obs; ltrend(ik)=ik; end
        for ik=1:obs; exosq(ik)=ik^2; end
        exo = [ones(obs,1) ltrend'];
        if sensi == 2
            exo = [ones(obs,1) ltrend' exosq'];
        end
        
        exo = exo(beg_samp:end_samp,:);
        
        samplelength=length(xdata(beg_samp:end_samp,:));

        %% estimate impulse responses
        [impzout,erzout,azeroout,a0betazout,var_explode]=...
            estimateIRF(xdata(beg_samp:end_samp,:),nlagq,nptsvar,teststatio,exo);       

        irf1(:,:,sensi,jrun) =  impzout(:,:,1);

        %% bootstrap IRF
        simdout=simdata(xdata(beg_samp:end_samp,:),...
        nlagq,nptsvar,ndraws,impzout,erzout,azeroout,a0betazout,exo);
        simimpzout=[];
        for zx=1:ndraws;
            [simimpzout,erzout,azeroout,a0betazout,var_explode] = estimateIRF(simdout(:,:,zx),...
            nlagq,nptsvar,teststatio,exo);                   
            irf1_mc(zx,:,:) = squeeze(simimpzout(:,:,1));
        end
    
        irfse1(:,:,sensi,jrun) = squeeze(std(irf1_mc(:,:,:),0,1));
    end
end


    %% plot
 for jrun=1:3;
 if jrun==1;
            namshort = char('g','y','c','lr','rx','dy','inf');       
            yshare = [gy 1 cy 1 1 1 1];
        elseif jrun ==2;
            namshort = char('g','y','i','lr','rx','dy','inf');       
            yshare = [gy 1 iy 1 1 1 1];
        elseif jrun ==3;
            namshort = char('g','y','nx','lr','rx','dy','inf');       
            yshare = [gy 1 1 1 1 1 1];
 end

     
        figure
        time = [0:nptsvar-1];
        for zx=1:nvar;
            subplot(3,3,zx)
            mycolors = [1 1 1; .85 .85 .85];            
            colormap(mycolors);

            rescale = max(irf1(1:10,1,1,jrun))*(yshare(1)/yshare(zx));
            point = irf1(:,zx,1,jrun)./rescale;
            sepoint = irfse1(:,zx,1,jrun)./rescale;
            ub = point+confidence*sepoint;
            lb = point-confidence*sepoint;  

            hold on;   
            gx=[lb ub-lb];
            area(time,gx,min(lb),'EdgeColor','none');
            plot(time,zeros(nptsvar,1),'k--','LineWidth',1);
         
            sensi = 2;
            rescale = max(irf1(1:10,1,sensi,jrun))*(yshare(1)/yshare(zx));
            point = irf1(:,zx,sensi,jrun)./rescale;
            plot(time,point,'r--','LineWidth',1);

            sensi = 3;
            rescale = max(irf1(1:10,1,sensi,jrun))*(yshare(1)/yshare(zx));
            point = irf1(:,zx,sensi,jrun)./rescale;
            plot(time,point,'k-.','LineWidth',1);

            sensi = 4;
            rescale = max(irf1(1:10,1,sensi,jrun))*(yshare(1)/yshare(zx));
            point = irf1(:,zx,sensi,jrun)./rescale;
            plot(time,point,'b-','LineWidth',1,'MarkerSize',2);

            
            axis([-.2 nptsvar-.8  min([1.5*min(lb) 0]) 1.5*max(ub)]);box on;       
            box on;
            title(deblank(namshort(zx,:)))        
        end
    end
    
if saveirf

    disp('writing figures to files...')
    for jrun=1:3;
        if jrun==1;
            namshort = char('g','y','c','lr','rx','dy','inf');       
            yshare = [gy 1 cy 1 1 1 1];
        elseif jrun ==2;
            namshort = char('g','y','i','lr','rx','dy','inf');       
            yshare = [gy 1 iy 1 1 1 1];
        elseif jrun ==3;
            namshort = char('g','y','nx','lr','rx','dy','inf');       
            yshare = [gy 1 1 1 1 1 1];
 end

    
    
    for zx=1:nvar;            
        figure
        set(gcf,'Visible','off')
        mycolors = [1 1 1; .85 .85 .85];            
        colormap(mycolors);

            rescale = max(irf1(1:10,1,1,jrun))*(yshare(1)/yshare(zx));
            point = irf1(:,zx,1,jrun)./rescale;
            sepoint = irfse1(:,zx,1,jrun)./rescale;
        ub = point+confidence*sepoint;
        lb = point-confidence*sepoint;  

        hold on;   
        gx=[lb ub-lb];
        area(time,gx,min(lb),'EdgeColor','none');
        plot(time,zeros(nptsvar,1),'k--','LineWidth',1);
        
           sensi = 2;
            rescale = max(irf1(1:10,1,sensi,jrun))*(yshare(1)/yshare(zx));
            point = irf1(:,zx,sensi,jrun)./rescale;
            plot(time,point,'r--','LineWidth',4);

            sensi = 3;
            rescale = max(irf1(1:10,1,sensi,jrun))*(yshare(1)/yshare(zx));
            point = irf1(:,zx,sensi,jrun)./rescale;
            plot(time,point,'k-.','LineWidth',4);

            sensi = 4;
            rescale = max(irf1(1:10,1,sensi,jrun))*(yshare(1)/yshare(zx));
            point = irf1(:,zx,sensi,jrun)./rescale;
            plot(time,point,'b-','LineWidth',4,'MarkerSize',2);
        
        
        axis([-.2 nptsvar-.8  min([1.5*min(lb) 0]) 1.5*max(ub)]);box on;       
        set(gca,'FontSize',25)
        box on;
        saveas(gcf,['..\figures\sensi_' num2str(jrun) '_' deblank(namshort(zx,:)) '.eps'],'psc2')
    end

end
end
    